HW_3

Data Import

data <- readRDS("very_low_birthweight.RDS")

Задание 1 1

Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего задания). Это данные о 671 младенце с очень низкой массой тела (<1600 грамм), собранные в Duke University Medical Center доктором Майклом О’Ши c 1981 по 1987 г.
Описание переменных см. здесь. Переменными исхода являются колонки ‘dead’, а также время от рождения до смерти или выписки (выводятся из ‘birth’ и ‘exit’ 7 пациентов были выписаны до рождения). Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.

data %>% 
  head()

##    birth   exit hospstay    lowph pltct  race  bwt gest        inout twn lol
## 1 81.511 81.604       34       NA   100 white 1250   35 born at Duke   0  NA
## 2 81.514 81.539        9 7.250000   244 white 1370   32 born at Duke   0  NA
## 3 81.552 81.552       -2 7.059998   114 black  620   23 born at Duke   0  NA
## 4 81.558 81.667       40 7.250000   182 black 1480   32 born at Duke   0  NA
## 5 81.593 81.599        2 6.969997    54 black  925   28 born at Duke   0  NA
## 6 81.602 81.771       62 7.189999    NA white  940   28 born at Duke   0  NA
##   magsulf meth toc  delivery apg1 vent pneumo pda cld      pvh      ivh    ipe
## 1      NA    0   0 abdominal    8    0      0   0   0     <NA>     <NA>   <NA>
## 2      NA    1   0 abdominal    7    0      0   0   0     <NA>     <NA>   <NA>
## 3      NA    0   1   vaginal    1    1      0   0  NA     <NA>     <NA>   <NA>
## 4      NA    1   0   vaginal    8    0      0   0   0     <NA>     <NA>   <NA>
## 5      NA    0   0 abdominal    5    1      1   0   0 definite definite   <NA>
## 6      NA    1   0 abdominal    8    1      0   0   0   absent   absent absent
##       year    sex dead
## 1 81.51196 female    0
## 2 81.51471 female    0
## 3 81.55304 female    1
## 4 81.55847   male    0
## 5 81.59406 female    1
## 6 81.60229 female    0

cols_to_select <- data %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "column", values_to = "NA_count") %>%
  filter(NA_count <= 100) %>%
  pull(column)

data_filtered = data %>% 
  select(all_of(cols_to_select)) %>% 
  na.omit()

data_filtered %>%
  write_csv("Low_BirthWeight_filtered_data.csv")

Задание 2

Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.

skim(data_filtered)

Data summary

Name	data_filtered
Number of rows	531
Number of columns	19
_______________________
Column type frequency:
factor	4
numeric	15
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
race	0	1	FALSE	4	bla: 303, whi: 211, nat: 13, ori: 4
inout	0	1	FALSE	2	bor: 448, tra: 83
delivery	0	1	FALSE	2	vag: 269, abd: 262
sex	0	1	FALSE	2	mal: 267, fem: 264

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
birth	0	1	84.63	1.54	81.51	83.43	84.77	85.83	87.48	▅▆▇▇▅
exit	0	1	84.76	1.55	81.05	83.56	84.87	85.99	87.72	▂▆▇▇▅
hospstay	0	1	47.04	63.50	-295.00	21.00	40.00	64.00	797.00	▁▇▁▁▁
lowph	0	1	7.22	0.13	6.53	7.13	7.22	7.32	7.55	▁▁▃▇▂
pltct	0	1	204.49	80.83	16.00	148.00	204.00	256.00	571.00	▂▇▅▁▁
bwt	0	1	1135.61	240.04	400.00	960.00	1160.00	1330.00	1500.00	▁▃▆▇▇
gest	0	1	29.25	2.21	23.00	28.00	29.00	31.00	36.00	▁▇▇▆▁
twn	0	1	0.21	0.40	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▂
apg1	0	1	5.02	2.65	0.00	2.00	5.00	7.00	9.00	▅▆▅▇▇
vent	0	1	0.54	0.50	0.00	0.00	1.00	1.00	1.00	▇▁▁▁▇
pneumo	0	1	0.18	0.38	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▂
pda	0	1	0.20	0.40	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▂
cld	0	1	0.26	0.44	0.00	0.00	0.00	1.00	1.00	▇▁▁▁▃
year	0	1	84.63	1.54	81.51	83.43	84.77	85.83	87.48	▅▆▇▇▅
dead	0	1	0.12	0.33	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▁

data_filtered %>% head()

##     birth   exit hospstay    lowph pltct  race  bwt gest        inout twn
## 2  81.514 81.539        9 7.250000   244 white 1370 32.0 born at Duke   0
## 4  81.558 81.667       40 7.250000   182 black 1480 32.0 born at Duke   0
## 5  81.593 81.599        2 6.969997    54 black  925 28.0 born at Duke   0
## 7  81.610 81.697       32 7.320000   282 black 1255 29.5 born at Duke   0
## 10 81.624 81.700       28 7.160000   153 black 1350 34.0 born at Duke   0
## 11 81.626 81.730       38 7.039997   229 white 1310 32.0 born at Duke   0
##     delivery apg1 vent pneumo pda cld     year    sex dead
## 2  abdominal    7    0      0   0   0 81.51471 female    0
## 4    vaginal    8    0      0   0   0 81.55847   male    0
## 5  abdominal    5    1      1   0   0 81.59406 female    1
## 7    vaginal    9    0      0   0   0 81.61053 female    0
## 10 abdominal    4    0      0   0   0 81.62421 female    0
## 11   vaginal    6    1      0   0   0 81.62695   male    0

data_filtered <- data_filtered %>% 
  mutate(across(where(is.factor), as.factor),
         across(where(~ is.numeric(.) && all(. %in% c(0, 1))), as.factor))

data_filtered %>% 
  select(where(is.numeric), inout) %>% 
  pivot_longer(cols = where(is.numeric) & !inout, names_to = "column", values_to = "value") %>% 
  ggplot(aes(x = value, fill = inout)) +
  geom_histogram(bins = 12, color = "black") +
  facet_wrap(~column, scales = "free_x") + 
  theme_minimal()

skim(data_filtered)

Data summary

Name	data_filtered
Number of rows	531
Number of columns	19
_______________________
Column type frequency:
factor	10
numeric	9
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
race	0	1	FALSE	4	bla: 303, whi: 211, nat: 13, ori: 4
inout	0	1	FALSE	2	bor: 448, tra: 83
twn	0	1	FALSE	2	0: 422, 1: 109
delivery	0	1	FALSE	2	vag: 269, abd: 262
vent	0	1	FALSE	2	1: 288, 0: 243
pneumo	0	1	FALSE	2	0: 438, 1: 93
pda	0	1	FALSE	2	0: 425, 1: 106
cld	0	1	FALSE	2	0: 393, 1: 138
sex	0	1	FALSE	2	mal: 267, fem: 264
dead	0	1	FALSE	2	0: 467, 1: 64

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
birth	0	1	84.63	1.54	81.51	83.43	84.77	85.83	87.48	▅▆▇▇▅
exit	0	1	84.76	1.55	81.05	83.56	84.87	85.99	87.72	▂▆▇▇▅
hospstay	0	1	47.04	63.50	-295.00	21.00	40.00	64.00	797.00	▁▇▁▁▁
lowph	0	1	7.22	0.13	6.53	7.13	7.22	7.32	7.55	▁▁▃▇▂
pltct	0	1	204.49	80.83	16.00	148.00	204.00	256.00	571.00	▂▇▅▁▁
bwt	0	1	1135.61	240.04	400.00	960.00	1160.00	1330.00	1500.00	▁▃▆▇▇
gest	0	1	29.25	2.21	23.00	28.00	29.00	31.00	36.00	▁▇▇▆▁
apg1	0	1	5.02	2.65	0.00	2.00	5.00	7.00	9.00	▅▆▅▇▇
year	0	1	84.63	1.54	81.51	83.43	84.77	85.83	87.48	▅▆▇▇▅

Удаление аутлайеров

Видим, что точно есть что-то странное в колонке hospstay - отрицательные значения. Посмотрим на них:

data_filtered %>% filter(hospstay < 0) %>% select(birth,exit, hospstay, dead) %>% mutate(hosp_calc = exit-birth)

##     birth   exit hospstay dead hosp_calc
## 1  81.826 81.826       -2    1     0.000
## 2  81.834 81.046     -288    0    -0.788
## 3  81.966 81.169     -291    0    -0.797
## 4  82.248 82.248       -2    1     0.000
## 5  83.116 83.116       -2    1     0.000
## 6  83.592 83.592       -2    1     0.000
## 7  83.732 83.732       -2    0     0.000
## 8  83.775 83.118     -240    1    -0.657
## 9  85.851 85.049     -293    0    -0.802
## 10 86.883 86.075     -295    0    -0.808
## 11 85.886 85.300     -214    0    -0.586

Вот -200 тут совсем странно выглядит, но не совсем понятно, почему у некоторых детей -2 стоит в значениях. МНе кажется, что это ошибка занесения и что на самом деле, тут должны быть 0 (как для большинства детей).

Еще видим, что у есть странные выбросы по pltct и lowph - низкие значения. Для их фильтрации возьмём 1.5 IQR метод (Значения должны быть между Q1-1.5IQR, Q3+1.5IQR).

UPD: потом нашёлся аутлайер по hospstay (>700 days) - его тоже пофильтруем по. Возьмём <400 days.

Итого:

data_filtered_outliers <- data_filtered %>% 
  filter(hospstay >= -3 & hospstay < 400,
         pltct < quantile(pltct, 0.75, na.rm = TRUE) + 1.5 * IQR(pltct, na.rm = TRUE) &
         pltct > quantile(pltct, 0.25, na.rm = TRUE) - 1.5 * IQR(pltct, na.rm = TRUE),
         lowph < quantile(lowph, 0.75, na.rm = TRUE) + 1.5 * IQR(lowph, na.rm = TRUE) &
         lowph > quantile(lowph, 0.25, na.rm = TRUE) - 1.5 * IQR(lowph, na.rm = TRUE)
          )

Все отфильтровали, было 531 наблюдение - стало 515.

data_filtered_outliers %>% 
  select(where(is.numeric), inout) %>% 
  pivot_longer(cols = where(is.numeric) & !inout, names_to = "column", values_to = "value") %>% 
  ggplot(aes(x = value, fill = inout)) +
  geom_histogram(bins = 12, color = "black") +
  facet_wrap(~column, scales = "free_x") + 
  theme_minimal()

# Задание 3

Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

library(rstatix)
data_filtered_outliers %>% 
  group_by(inout) %>% 
  select(inout, lowph) %>% 
  dfSummary()

## Data Frame Summary  
## data_filtered_outliers  
## Group: inout = born at Duke  
## Dimensions: 433 x 2  
## Duplicates: 376  
## 
## --------------------------------------------------------------------------------------------------------
## No   Variable    Stats / Values          Freqs (% of Valid)   Graph                 Valid      Missing  
## ---- ----------- ----------------------- -------------------- --------------------- ---------- ---------
## 2    lowph       Mean (sd) : 7.2 (0.1)   57 distinct values           :   :         433        0        
##      [numeric]   min < med < max:                                     : . :         (100.0%)   (0.0%)   
##                  6.9 < 7.2 < 7.5                                    . : : : .                           
##                  IQR (CV) : 0.2 (0)                                 : : : : :                           
##                                                               . . : : : : : : . .                       
## --------------------------------------------------------------------------------------------------------
## 
## Group: inout = transported  
## Dimensions: 80 x 2  
## Duplicates: 43  
## 
## --------------------------------------------------------------------------------------------------------
## No   Variable    Stats / Values          Freqs (% of Valid)   Graph                 Valid      Missing  
## ---- ----------- ----------------------- -------------------- --------------------- ---------- ---------
## 2    lowph       Mean (sd) : 7.2 (0.1)   37 distinct values             : .         80         0        
##      [numeric]   min < med < max:                                       : :         (100.0%)   (0.0%)   
##                  6.9 < 7.2 < 7.4                                        : :                             
##                  IQR (CV) : 0.2 (0)                             . :   : : : : : :                       
##                                                               . : : . : : : : : :                       
## --------------------------------------------------------------------------------------------------------

data_filtered_outliers %>%
  t_test(lowph ~ inout) %>%
  add_significance()

## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "born at Duke" - "transported", or divided in the order "born at Duke" /
## "transported" for ratio-based statistics. To specify this order yourself,
## supply `order = c("born at Duke", "transported")`.

## # A tibble: 1 × 7
##   statistic  t_df     p_value alternative estimate lower_ci upper_ci
##       <dbl> <dbl>       <dbl> <chr>          <dbl>    <dbl>    <dbl>
## 1      5.69  108. 0.000000108 two.sided     0.0847   0.0552    0.114

data_filtered_outliers %>% select(inout, lowph) %>%
  pivot_longer(cols = lowph, names_to='lowph', values_to='value') %>%
  ggplot(aes(x = lowph, y = value, fill=inout)) +
  geom_boxplot() +
  theme_minimal()

Анализ показал, что уровень рН у пациентов, транспортированных из других учреждений, значительно ниже, чем у новорожденных в госпитале Duke (MD = 0.084, 95% CI: [0.055; 0.114], p < 0.0001). Поскольку, предположительно, снижение уровня рН ассоциируется с повышенной смертностью, можно сделать вывод, что транспортировка пациентов может быть потенциальным фактором риска, способствующим увеличению вероятности летального исхода.

Задание 4

Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

corr_data <- data_filtered_outliers %>% 
  select(where(is.numeric), -c(birth, year, exit))

сorr_map = cor(corr_data)
corrplot(сorr_map, method = "color", type = "lower", 
         addCoef.col = "grey30", diag = FALSE,
         cl.pos = "b", tl.col = "grey10",
         col = COL2('RdBu', 10))

library(corrr)

## Warning: package 'corrr' was built under R version 4.4.2

## 
## Attaching package: 'corrr'

## The following object is masked from 'package:skimr':
## 
##     focus

сorr_map %>% 
  network_plot(min_cor = .0)

Время беременности коррелирует с весом плода.

Задание 5

Постройте иерархическую кластеризацию на этом датафрейме.

data_filtered_outliers <- data_filtered_outliers %>% 
  mutate(ID = c(1:nrow(data_filtered_outliers)))

rownames(data_filtered_outliers) <- data_filtered_outliers$ID

res <- data_filtered_outliers %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  scale() %>% 
  dist(method = "euclidean") %>% 
  hclust(method = "ward.D2") 
res %>% 
  fviz_dend(cex = 0.2, 
            k = 4, 
            k_colors = "jco")

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Нашёлся ещё какой-то аутлайер - попал в отдельный кластер. UPD там по hospstay был. Пофильтровал выше и поменял графики. Получили 4 кластера.

Задание 6

Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.

data_filtered_outliers %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  scale() %>% 
  pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50), 
                     fontsize_row = 0.00001, fontsize_col = 10, angle_col = 0)

data_filtered_outliers %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>%
  scale() %>% 
  pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50), 
                     fontsize_row = 0.00001, fontsize_col = 8, angle_col = 0, kmeans_k=4)

Пользуясь тем, что в иерархической кластеризации мы остановились на 4 кластерах. Тут тоже сфокусируемся на 4 кластерах. Видим, что есть кластер (самый низкий) с очень высоким hospstay и низким bwt => что логично, дети с низким весом требуют ухода и присмотра врачей.

Первый кластер (самый высокий) с самым низким относительный hospstay, высоким bwt, lowph, pltct => дети, которые были +- здоровые, поэтому быстро были выписаны.

Два кластера посредеине со средним hospstay, но различающиеся по средним уровням клинических метрик для детей. Вот эти кластера, по-хорошему, надо бы ещё поделить и посмотреть, что в них внутри.

Задание 7

Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?

В нашем случае: во-первых, удалим не нужные колонки, которые являются датой, которая нам ничего не даст в рамках EDA; во-вторых, шкалируем, т.к. скейлы колонок разные!

pca_data = data_filtered_outliers %>% 
  select(where(is.numeric), -c(ID, birth, year, exit)) %>% 
  scale() %>%
  prcomp(scale=TRUE)

fviz_eig(pca_data, 
         addlabels = T, 
         ylim = c(0, 50))

pca_scores <- as.data.frame(pca_data$x)

ggplot(pca_scores, aes(x = PC1, y = PC2)) +
  geom_point(color = "steelblue", size = 2) +
  theme_minimal() +
  labs(
    title = "PCA Plot",
    x = "Principal Component 1",
    y = "Principal Component 2"
  )

Видим, что сильных аутлайеров у нас нет. Первые две компоненты объясняют 60% дисперсии в данных, что достаточно много.

Задание 8

Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.

ggbiplot(pca_data, choices=1:2,
         scale=1, alpha = 0.5, groups = as.factor(data_filtered_outliers$dead)) + 
  theme_bw()

У нас много корреляций между hospstay и остальными числовыми колонками, поэтому первая компоненты такая сильная.

Задание 9

Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.

pca_scores <- as.data.frame(pca_data$x) %>%
  mutate(ID = data_filtered_outliers$ID,
         group = as.factor(data_filtered_outliers$dead))

library(plotly)

## Warning: package 'plotly' was built under R version 4.4.2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

plot_ly(
  data = pca_scores,
  x = ~PC1,
  y = ~PC2,
  type = 'scatter',
  mode = 'markers',
  color = ~group,
  text = ~paste('ID:', ID),  # Текст для отображения при наведении
  hoverinfo = 'text'
) %>%
  layout(
    title = 'PCA Plot (PC1 vs PC2)',
    xaxis = list(title = 'Principal Component 1'),
    yaxis = list(title = 'Principal Component 2'),
    legend = list(title = list(text = 'Group'))
  )

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Задание 10

Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?

Используя свой опыт из онкологии могут сказать следущее: события dead = 1 (смерть), dead = 0 (цензурирование), могли произойти по многим причинам: ребенку действительно было плохо, халатность врачей, случайная трагедия и т.д. Поэтому для оценки выживаемости колонку исхода всегда используют вместе с самой выживаемостью (в днях, например).

Задание 11

Приведите ваши данные к размерности в две колонки через UMAP. Сравните результаты отображения точек между алгоритмами PCA и UMAP

library(embed)

umap_data = data_filtered_outliers %>% select(where(is.numeric), -c(ID, birth, year, exit))
umap_prep <- recipe(~., data = umap_data) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>%
  step_umap(all_predictors()) %>%
  prep() %>%
  juice()

umap_scores <- as.data.frame(umap_prep) %>%
  mutate(group = as.factor(data_filtered_outliers$dead))

umap_scores %>%
  ggplot(aes(x=UMAP1, y=UMAP2, color=group)) + #  # можно добавить раскраску 
  geom_point(alpha = 0.7, size = 2) +
  labs(color = NULL) +
  theme_minimal()

Задание 12

Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Измените основные параметры UMAP (n_neighbors и min_dist) и проанализируйте, как это влияет на результаты.

umap_data = data_filtered_outliers %>% select(where(is.numeric), -c(ID, birth, year, exit))
umap_prep <- recipe(~., data = umap_data) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>%
  step_umap(all_predictors(), min_dist=0.001, neighbors=5) %>%
  prep() %>%
  juice()

umap_scores <- as.data.frame(umap_prep) %>%
  mutate(group = as.factor(data_filtered_outliers$dead))

umap_scores %>%
  ggplot(aes(x=UMAP1, y=UMAP2, color=group)) + #  # можно добавить раскраску 
  geom_point(alpha = 0.7, size = 2) +
  labs(color = NULL) +
  theme_minimal()

umap_data = data_filtered_outliers %>% select(where(is.numeric), -c(ID, birth, year, exit))
umap_prep <- recipe(~., data = umap_data) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>%
  step_umap(all_predictors(), min_dist=0.5, neighbors=15) %>%
  prep() %>%
  juice()

umap_scores <- as.data.frame(umap_prep) %>%
  mutate(group = as.factor(data_filtered_outliers$dead))

umap_scores %>%
  ggplot(aes(x=UMAP1, y=UMAP2, color=group)) + #  # можно добавить раскраску 
  geom_point(alpha = 0.7, size = 2) +
  labs(color = NULL) +
  theme_minimal()

По мои оценкам: neighbors сильнее влияет на график UMAP, чем min_dist. Neighbors влияет на “кластеризацию точек”, а мин-дист больше на их разброс.